R Notebook

Dataset ejemplo: lung (paquete survival)

#| echo: true
library(survival)
data(cancer, package="survival")
table(lung$ph.ecog)

  0   1   2   3 
 63 113  50   1 

lung <- lung %>% 
  na.omit() %>%
  mutate(status2 = ifelse(status == 2, 1, 0)) %>%
  filter(ph.ecog<3) %>%
  mutate(surv_obj = Surv(time, status2),
         sex = factor(sex, 
                      levels = c(1, 2), 
                      labels = c("Hombre", "Mujer")),
         ph.ecog2 = factor(ph.ecog, 
                           levels = c(0,1,2), 
                           labels = c("asymptomatic", 
                                   "symptomatic but completely ambulatory",
                                   "in bed <50% of the day")))
  

head(lung)

summary(lung)

      inst            time            status           age            sex         ph.ecog          ph.karno     
 Min.   : 1.00   Min.   :   5.0   Min.   :1.000   Min.   :39.00   Hombre:102   Min.   :0.0000   Min.   : 50.00  
 1st Qu.: 3.00   1st Qu.: 175.2   1st Qu.:1.000   1st Qu.:57.00   Mujer : 64   1st Qu.:0.0000   1st Qu.: 70.00  
 Median :11.00   Median : 268.5   Median :2.000   Median :63.50                Median :1.0000   Median : 80.00  
 Mean   :10.69   Mean   : 311.1   Mean   :1.717   Mean   :62.52                Mean   :0.9458   Mean   : 82.17  
 3rd Qu.:15.00   3rd Qu.: 422.8   3rd Qu.:2.000   3rd Qu.:70.00                3rd Qu.:1.0000   3rd Qu.: 90.00  
 Max.   :32.00   Max.   :1022.0   Max.   :2.000   Max.   :82.00                Max.   :2.0000   Max.   :100.00  
   pat.karno         meal.cal         wt.loss           status2          surv_obj.time       surv_obj.status  
 Min.   : 30.00   Min.   :  96.0   Min.   :-24.000   Min.   :0.0000   Min.   :   5.0000    Min.   :0.0000000  
 1st Qu.: 70.00   1st Qu.: 616.0   1st Qu.:  0.000   1st Qu.:0.0000   1st Qu.: 175.2500    1st Qu.:0.0000000  
 Median : 80.00   Median : 975.0   Median :  6.500   Median :1.0000   Median : 268.5000    Median :1.0000000  
 Mean   : 79.64   Mean   : 928.2   Mean   :  9.657   Mean   :0.7169   Mean   : 311.0904    Mean   :0.7168675  
 3rd Qu.: 90.00   3rd Qu.:1168.8   3rd Qu.: 15.000   3rd Qu.:1.0000   3rd Qu.: 422.7500    3rd Qu.:1.0000000  
 Max.   :100.00   Max.   :2600.0   Max.   : 68.000   Max.   :1.0000   Max.   :1022.0000    Max.   :1.0000000  
                                  ph.ecog2 
 asymptomatic                         :47  
 symptomatic but completely ambulatory:81  
 in bed <50% of the day               :38  
                                           
                                           
                                           

#| echo: true
print(lung$surv_obj)

  [1]  455   210  1022+  310   361   218   166   170   567   613   707    61   301    81   371   520   574   390    12 
 [20]  473    26   107    53   814   965+   93   731   460   153   433   583    95   303   519   643   765    53   246 
 [39]  689     5   687   345   444   223    60   163    65   821+  428   230   840+  305    11   226   426   705   363 
 [58]  176   791    95   196+  167   806+  284   641   147   740+  163   655    88   245    30   477   559+  450   156 
 [77]  529+  429   351    15   181   283    13   212   524   288   363   199   550    54   558   207    92    60   551+
 [96]  293   353   267   511+  457   337   201   404+  222    62   458+  353   163    31   229   156   291   179   376+
[115]  384+  268   292+  142   413+  266+  320   181   285   301+  348   197   382+  303+  296+  180   145   269+  300+
[134]  284+  292+  332+  285   259+  110   286   270   225+  269   225+  243+  276+  135    79    59   240+  202+  235+
[153]  239   252+  221+  185+  222+  183   211+  175+  197+  203+  191+  105+  174+  177+

Evaluación inicial: ¿sex viola PH?

mod1 <- coxph(surv_obj ~ sex + age, data = lung)
summary(mod1)

Call:
coxph(formula = surv_obj ~ sex + age, data = lung)

  n= 166, number of events= 119 

             coef exp(coef) se(coef)      z Pr(>|z|)  
sexMujer -0.43909   0.64462  0.19791 -2.219   0.0265 *
age       0.01679   1.01694  0.01086  1.546   0.1221  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

         exp(coef) exp(-coef) lower .95 upper .95
sexMujer    0.6446     1.5513    0.4374    0.9501
age         1.0169     0.9833    0.9955    1.0388

Concordance= 0.597  (se = 0.032 )
Likelihood ratio test= 8.43  on 2 df,   p=0.01
Wald test            = 8.05  on 2 df,   p=0.02
Score (logrank) test = 8.17  on 2 df,   p=0.02

Evaluemos el supuesto de Proporcionalidad

zph1 <- cox.zph(mod1)
zph1

       chisq df    p
sex    0.912  1 0.34
age    0.557  1 0.46
GLOBAL 1.296  2 0.52

plot(zph1)

• Si el p-valor en la prueba cox.zph() para sex es significativo o el gráfico tiene tendencia → violación de PH.

Estimador de Kaplan-Meier por sexo

#| echo: true
fit_km <- survfit(surv_obj ~ sex, data = lung)

ggsurvplot(fit_km,
           data = lung,
           conf.int = TRUE,
           xlab = "Días", ylab = "Supervivencia estimada",
           title = "Curvas Kaplan-Meier por sexo")

ggsurvplot(fit_km,
           data = lung,
           fun = "cloglog",
           conf.int = TRUE,
           xlab = "Días", ylab = "Supervivencia estimada",
           title = "Curvas Kaplan-Meier por sexo")

Estratificación por sex

#| echo: true
mod_strat <- coxph(surv_obj ~ age + strata(sex), data = lung)
summary(mod_strat)

Call:
coxph(formula = surv_obj ~ age + strata(sex), data = lung)

  n= 166, number of events= 119 

       coef exp(coef) se(coef)     z Pr(>|z|)
age 0.01600   1.01613  0.01085 1.475     0.14

    exp(coef) exp(-coef) lower .95 upper .95
age     1.016     0.9841    0.9948     1.038

Concordance= 0.558  (se = 0.031 )
Likelihood ratio test= 2.23  on 1 df,   p=0.1
Wald test            = 2.18  on 1 df,   p=0.1
Score (logrank) test = 2.18  on 1 df,   p=0.1

zph_strat <- cox.zph(mod_strat)
zph_strat

       chisq df    p
age    0.503  1 0.48
GLOBAL 0.503  1 0.48

plot(zph_strat)

• La función de riesgo base es diferente para hombres y mujeres.
• La covariable age tiene el mismo efecto en ambos estratos.

Comparación de modelos: sin estratificación vs con estratificación

#| echo: true
# Modelo sin estratificación
AIC(mod1)

[1] 1001.58

# Modelo con estratificación
AIC(mod_strat)

[1] 854.7006

# Ver ambos resúmenes
summary(mod1)

Call:
coxph(formula = surv_obj ~ sex + age, data = lung)

  n= 166, number of events= 119 

             coef exp(coef) se(coef)      z Pr(>|z|)  
sexMujer -0.43909   0.64462  0.19791 -2.219   0.0265 *
age       0.01679   1.01694  0.01086  1.546   0.1221  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

         exp(coef) exp(-coef) lower .95 upper .95
sexMujer    0.6446     1.5513    0.4374    0.9501
age         1.0169     0.9833    0.9955    1.0388

Concordance= 0.597  (se = 0.032 )
Likelihood ratio test= 8.43  on 2 df,   p=0.01
Wald test            = 8.05  on 2 df,   p=0.02
Score (logrank) test = 8.17  on 2 df,   p=0.02

#| echo: true
summary(mod_strat)

Call:
coxph(formula = surv_obj ~ age + strata(sex), data = lung)

  n= 166, number of events= 119 

       coef exp(coef) se(coef)     z Pr(>|z|)
age 0.01600   1.01613  0.01085 1.475     0.14

    exp(coef) exp(-coef) lower .95 upper .95
age     1.016     0.9841    0.9948     1.038

Concordance= 0.558  (se = 0.031 )
Likelihood ratio test= 2.23  on 1 df,   p=0.1
Wald test            = 2.18  on 1 df,   p=0.1
Score (logrank) test = 2.18  on 1 df,   p=0.1

anova(mod1, mod_strat, test = "LRT")

Analysis of Deviance Table
 Cox model: response is  surv_obj
 Model 1: ~ sex + age
 Model 2: ~ age + strata(sex)
   loglik  Chisq Df Pr(>|Chi|)    
1 -498.79                         
2 -426.35 144.88  1  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

• Se puede comparar el AIC para ver si mejora el ajuste.
• Nota: el coeficiente de sex desaparece en el modelo estratificado.

Comparación visual de curvas ajustadas

#| echo: true
library(survminer)
fit_strat <- survfit(mod_strat)

ggsurvplot(fit_strat, data = lung,
           legend.title = "Sexo (estratos)",
           xlab = "Días", ylab = "Supervivencia ajustada")

¿Qué sigue?

#| echo: true
mod_strat2 <- coxph(surv_obj ~ age + sex + ph.karno + pat.karno + meal.cal + 
                      ph.ecog2 + wt.loss, data = lung)
summary(mod_strat2)

Call:
coxph(formula = surv_obj ~ age + sex + ph.karno + pat.karno + 
    meal.cal + ph.ecog2 + wt.loss, data = lung)

  n= 166, number of events= 119 

                                                    coef  exp(coef)   se(coef)      z Pr(>|z|)   
age                                            9.968e-03  1.010e+00  1.173e-02  0.849  0.39563   
sexMujer                                      -5.489e-01  5.776e-01  2.023e-01 -2.713  0.00667 **
ph.karno                                       2.213e-02  1.022e+00  1.133e-02  1.954  0.05072 . 
pat.karno                                     -1.146e-02  9.886e-01  8.579e-03 -1.336  0.18160   
meal.cal                                       2.708e-05  1.000e+00  2.609e-04  0.104  0.91733   
ph.ecog2symptomatic but completely ambulatory  6.474e-01  1.911e+00  2.822e-01  2.294  0.02178 * 
ph.ecog2in bed <50% of the day                 1.454e+00  4.282e+00  4.621e-01  3.147  0.00165 **
wt.loss                                       -1.430e-02  9.858e-01  7.830e-03 -1.827  0.06775 . 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                                              exp(coef) exp(-coef) lower .95 upper .95
age                                              1.0100     0.9901    0.9871    1.0335
sexMujer                                         0.5776     1.7314    0.3885    0.8587
ph.karno                                         1.0224     0.9781    0.9999    1.0453
pat.karno                                        0.9886     1.0115    0.9721    1.0054
meal.cal                                         1.0000     1.0000    0.9995    1.0005
ph.ecog2symptomatic but completely ambulatory    1.9106     0.5234    1.0989    3.3217
ph.ecog2in bed <50% of the day                   4.2822     0.2335    1.7311   10.5931
wt.loss                                          0.9858     1.0144    0.9708    1.0010

Concordance= 0.651  (se = 0.029 )
Likelihood ratio test= 26.72  on 8 df,   p=8e-04
Wald test            = 26.35  on 8 df,   p=9e-04
Score (logrank) test = 27.34  on 8 df,   p=6e-04

zph_strat2 <- cox.zph(mod_strat2)
zph_strat2

            chisq df     p
age        0.0541  1 0.816
sex        1.4949  1 0.221
ph.karno   5.5965  1 0.018
pat.karno  3.2368  1 0.072
meal.cal   6.6377  1 0.010
ph.ecog2   3.7785  2 0.151
wt.loss    0.0860  1 0.769
GLOBAL    14.3160  8 0.074

plot(zph_strat2)

modelo_completo <- coxph(surv_obj ~ age + strata(sex) + ph.karno + pat.karno + meal.cal + 
                      ph.ecog2 + wt.loss, data = lung)
birthwt.step <- stepAIC(modelo_completo, trace = TRUE)

Start:  AIC=848.42
surv_obj ~ age + strata(sex) + ph.karno + pat.karno + meal.cal + 
    ph.ecog2 + wt.loss

            Df    AIC
- meal.cal   1 846.43
- age        1 846.91
- pat.karno  1 848.36
<none>         848.42
- ph.karno   1 849.72
- wt.loss    1 850.36
- ph.ecog2   2 853.97

Step:  AIC=846.43
surv_obj ~ age + strata(sex) + ph.karno + pat.karno + ph.ecog2 + 
    wt.loss

            Df    AIC
- age        1 844.94
<none>         846.43
- pat.karno  1 846.43
- ph.karno   1 847.74
- wt.loss    1 848.37
- ph.ecog2   2 851.99

Step:  AIC=844.94
surv_obj ~ strata(sex) + ph.karno + pat.karno + ph.ecog2 + wt.loss

            Df    AIC
<none>         844.94
- pat.karno  1 845.00
- ph.karno   1 845.86
- wt.loss    1 847.13
- ph.ecog2   2 850.50

#| echo: true
mod_select <- coxph(surv_obj ~ strata(sex) + pat.karno + ph.karno + ph.ecog2 + wt.loss, data = lung)
summary(mod_select)

Call:
coxph(formula = surv_obj ~ strata(sex) + pat.karno + ph.karno + 
    ph.ecog2 + wt.loss, data = lung)

  n= 166, number of events= 119 

                                                   coef exp(coef)  se(coef)      z Pr(>|z|)   
pat.karno                                     -0.012249  0.987826  0.008474 -1.445  0.14833   
ph.karno                                       0.018725  1.018901  0.011237  1.666  0.09565 . 
ph.ecog2symptomatic but completely ambulatory  0.583929  1.793070  0.282687  2.066  0.03886 * 
ph.ecog2in bed <50% of the day                 1.419825  4.136396  0.464084  3.059  0.00222 **
wt.loss                                       -0.015643  0.984479  0.007915 -1.976  0.04813 * 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                                              exp(coef) exp(-coef) lower .95 upper .95
pat.karno                                        0.9878     1.0123    0.9716    1.0044
ph.karno                                         1.0189     0.9814    0.9967    1.0416
ph.ecog2symptomatic but completely ambulatory    1.7931     0.5577    1.0303    3.1205
ph.ecog2in bed <50% of the day                   4.1364     0.2418    1.6657   10.2719
wt.loss                                          0.9845     1.0158    0.9693    0.9999

Concordance= 0.636  (se = 0.03 )
Likelihood ratio test= 19.99  on 5 df,   p=0.001
Wald test            = 20.75  on 5 df,   p=9e-04
Score (logrank) test = 21.54  on 5 df,   p=6e-04

zph_select <- cox.zph(mod_select)
zph_select

           chisq df     p
pat.karno 2.4778  1 0.115
ph.karno  5.1171  1 0.024
ph.ecog2  3.7148  2 0.156
wt.loss   0.0602  1 0.806
GLOBAL    6.5404  5 0.257

plot(zph_select)